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New generation in vitro high-throughput screening (HTS) assays 
for the assessment of engineered nanomaterials provide an opportu- 
nity to learn how these particles interact at the cellular level, par- 
ticularly in relation to injury pathways. These types of assays are 
often characterized by small sample sizes, high measurement error 
and high dimensionality, as multiple cytotoxicity outcomes are mea- 
sured across an array of doses and durations of exposure. In this 
paper we propose a probability model for the toxicity profiling of en- 
gineered nanomaterials. A hierarchical structure is used to account 
for the multivariate nature of the data by modeling dependence be- 
tween outcomes and thereby combining information across cytotox- 
icity pathways. In this framework we are able to provide a flexible 
surface-response model that provides inference and generalizations of 
various classical risk assessment parameters. We discuss applications 
of this model to data on eight nanoparticles evaluated in relation to 
four cytotoxicity parameters. 

1. Introduction. Nanotechnology is rapidly growing and currently used 
in various industries such as food, agriculture, electronics, textiles and health 
care. The widespread use of engineered nanomaterials (ENM) in over 800 
consumer products increases the likelihood that these materials will come 
into contact with humans and the environment [Maynard et al. (2006), 
Kahru and Dubourguier (2009)]. Many biological processes take place at 
the nanoscale level, and the introduction of ENMs into living organisms 
could lead to interference in the molecular and cellular processes that are 
critical to life [Nel et al. (2009)]. This potential for human and environmen- 
tal hazard has spurred recent interest in early identification of potentially 
hazardous nanomaterials. Knowledge about the potential hazard of nanoma- 
terials is still lacking and a lot of study is required to understand how ENM 
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Fig. 1. Fluorescence images and heat map of raw data. On the left are fluorescence 
images of RAW cells treated with various nanomaterials (quantum dot, platinum and a 
negative control) and dyed with compatible dye combinations including MitoSox, JC1. PI 
and Fluo-4- The subsequence fluorescence read-out, measured at varying wavelengths, pro- 
vides a measure of the number of cells positive for the response. On the left is a heat map 
of the raw data for each particle and outcome. Colder colors indicate a smaller percentage 
of cells positive for the response and warmer colors indicate a larger percentage of cells 
positive for the response. 

properties such as size, shape, agglomeration state, solubility and surface 
properties could lead to hazard generation at the nano-bio interface [Stern 
and McNeil (2008), Nel et al. (2006)]. 

Current research in nano-toxicology includes new generation high-through- 
put screening (HTS) assays, which enable the simultaneous observation of 
multiple cellular injury pathways across an array of doses and times of ex- 
posure. In this article, for example, we analyze data on eight metal and 
metal oxide nanoparticles, monitored in relation to four cellular injury re- 
sponses, derived from the hierarchical oxidative stress model of Nel et al. 
(2006) and Xia et al. (2006). All four outcomes are measured contempora- 
neously over a grid of ten doses and seven hours of exposure (see Figure 1). 
The four measured responses include mitochondrial superoxide formation, 
loss of mitochondrial membrane potential, elevated intracellular calcium and 
membrane damage [George et al. (2009)]. For increasing dosage and dura- 
tion of exposure, we observe typical dose-response kinetics, with outcomes 
possibly depending on one another. 

These assays provide an opportunity to help define biological relation- 
ships and may suggest which nanoparticles are likely to have an in vivo 
effect. While HTS assays cannot replace traditional animal studies, they are 
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less costly, less labor intensive and can be used to explore the large number 
of potential nanomaterial variables that can influence human health hazards 
[Meng et al. (2010), Stanley et al. (2008), Maynard et al. (2006)]. The feasi- 
bility and utility of HTS assays have been illustrated in various fields such 
as functional genomics, with the use of microarray technology, as well as 
in pharmacology for the rapid screening of potential drug targets [Hoheisel 
(2006), White (2000)]. In toxicology, risk assessment involves the character- 
ization of hazard as well as the potential for exposure while accounting for 
all assumptions and uncertainties. The HTS framework provides a wealth 
of information about cellular injury pathways but proves a challenge for the 
classic risk assessment paradigm. In fact, there is still disagreement in the 
HTS setting on how to define and how appropriate are classical risk assess- 
ment parameters such as no observable adverse effect level (NOAEL), the 
lowest observable adverse effect level (LOAEL) and the dose that produces 
50% of the maximum response (EC50), among others. 

Parametric functions such as families of sigmoidal curves are frequently 
used to fit dose-response data. Some commonly used sigmoidal models in- 
clude log-logistic models, log-normal models and Weibull models [see Ritz 
(2010) for a recent review of these models]. The log-logistic functions are 
the most frequently used for modeling dose-response data in toxicology. The 
four parameter log-logistic model can be expressed as follows: 

d-c 

(!) f{x;b,c,d,h) = c+— — - — — — - — — — . 

1 + exp[b{log(x) - log(ft)}] 

Here h, the inflection point in the curve, provides a convenient risk assess- 
ment parameter, since it can be interpreted as the 50% effective or inhibitory 
dose (EC50, IC50) [Emmens (1940)]. Other special cases of this model in- 
clude the 3 parameter log-logistic model which leads to the famous Hill 
equation [Hill (1910)] and special cases of the Michaelis-Menton kinetics. 
Further extensions of these models include the five parameter log-logistic 
function, which provides a bit more flexibility by allowing the function to be 
asymmetric [Finney (1979)], and the Brain-Cousens model, which includes 
an extra parameter to account for a possible favorable response to a toxin at 
low concentrations [Calabrese and Baldwin (2003)]. In general, these models 
assume that the dose-response function is completely known apart from the 
few parameters to be estimated, usually by determining which values of the 
parameters result in the best fit to the dose-response function. 

Several other methods have been proposed to model nonlinear dose-response 
relationships relaxing strictly parametric assumptions. Ramsay (1988) pro- 
posed the use of monotone regression splines to model a dose-response func- 
tion. In this case, piecewise polynomials or splines can allow greater flexibil- 
ity while achieving monotonicity by imposing constraints on the estimated 
function. Li and Hunt (2004) proposed the use of linear B-splines with one 
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random interior knot to model a nonlinear dose-response curve. In this con- 
text, the random interior knot provides inference on the dose at which the 
toxin begins to take effect and thereby provides a useful parameter for risk 
assessment. Kong and Eubank (2006) suggested the use of functions that 
combine smoothing spline techniques and the nonnegativity properties of 
cubic B-splines to estimate the dose-response curve. The use of nonpar a- 
metric techniques to estimate dose-response curves often provides a more 
realistic representation of the data generating process. At the same time, 
however, some of these techniques make it more difficult to interpret the 
model in terms of classical risk assessment. 

Recent literature advocates the simultaneous use of multiple outcomes to 
assess risk. Regan and Catalano (1999) proposed a bivariate dose-response 
model that accounts for the dependence among outcomes of developmental 
toxicity using generalized estimating equations. Geys et al. (2001) proposed a 
similar model for risk assessment of developmental toxicity, but approached 
the problem using latent variables. Yu and Catalano (2005) suggested a 
model for quantitative risk assessment of bivariate continuous measures of 
neurotoxicity using percentile regression. These methods are often aimed at 
the analysis of one potentially toxic agent as it relates to adverse events or 
continuous outcomes observed in association with exposure over a range of 
doses. Their direct applicability to the general HTS setting described earlier 
is therefore limited. 

From a statistical perspective, cellular interrogation data based on high- 
throughput platforms can be characterized as multivariate dependent ob- 
servations. Each nanoparticle is indeed associated with a multiple set of 
cellular outcomes recorded both longitudinally, in relation to different expo- 
sure durations, and cross-sectionally, in relation to a dose escalation design. 
This particular design structure suggests that valid statistical inference must 
account for potentially complex patterns of dependence between different 
observations. A reasonable dependence scheme might, for example, assume 
data to be dependent within outcome and particle, as well as between out- 
comes for the same particle. 

In conjunction with considerations related to the joint sampling distribu- 
tion of these data structures, appropriate statistical treatment must account 
for nonlinearities in the mean response associated with dose and duration 
dynamics. While, in principle, one can choose to define a random response 
surface in a completely nonparametric fashion, it is important to maintain 
a certain degree of interpretability, especially in relation to standard haz- 
ard assessment quantities of interest to substantive scientists. In summary, 
perhaps reductively, the overall modeling challenge lies in the definition of a 
flexible and interpretable probabilistic representation for a family of depen- 
dent dose-response random surfaces. 

In this paper we propose a hierarchical dose-response model for the anal- 
ysis of HTS data from nanotoxicology. Our model builds on earlier work 
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[Hastie and Tibshirani (1986), Li and Hunt (2004)], expanding on them to 
account for the multivariate nature of the data and to address the estimation 
of a series of two-dimensional dose-response surfaces. We provide a flexible 
framework for modeling dose and duration response kinetics jointly, while 
providing inference on several risk assessment parameters of interest. We 
utilize a hierarchical structure to define dependence between outcomes and 
thereby borrow strength across injury pathways, providing the basis for a 
comprehensive risk assessment paradigm in HTS studies. We account for 
outlying observations via a T-distributed error model and describe how to 
carry out inference for the model parameters and their functions on the ba- 
sis of simulated draws from their posterior distribution. To our knowledge, 
we are the first to propose a principled statistical methodology for the joint 
analysis of this new generation of in vitro data. 

The remainder of the article is organized as follows. In Section 2 we intro- 
duce the proposed model. In Section 3 we discuss parameter estimation and 
associated inferential details. Section 4 employs the proposed model for the 
analysis of 8 metal oxide nanomaterials and describes inference for various 
risk assessment parameters of interest. We conclude with a critical discussion 
of the limitations and possible extensions of our method in Section 5. 

2. Model formulation. 

2.1. Model description. In this section we describe a dose-response model 
for a general HTS study, where we monitor a multivariate continuous out- 
come y, corresponding to J cytotoxicity parameters, in association with 
the exposure of a number of cells to / different ENMs. More precisely, 
let Uijk(d,t) denote a multivariate response corresponding to ENM i (i = 
1, . . . , I), cytotoxicity parameter j (j = 1, . . . , J) and replicate k (k = 1, . . . , K) 
at dose d G [0,D] and time t € [0,T]. In typical applications one observes 
y over a discrete set of doses d= (d\, . . . , d mi )' and exposure times t = 

(ii, ,t m2 )'. However, for clarity of exposition, we simplify our notation 

and without loss of generality refer to a general dose d € [0, D] and time 
t € [0,T]. We introduce the following 4-stage hierarchical model. 

Stage 1: Sampling model The observed response of particle i, cytotoxicity 
parameter j and replicate k is modeled as 



where Sijk(d,t) ~ iV(0, of./Tj). Here mij(d,t) denotes the response surface 
for particle i and outcome j. The proposed response surface describes dose 
and duration kinetics for all d £ [0, D] and t € [0, T] and is expected to exhibit 
a nonlinear dynamic over these domains. The distribution of y^ is modeled 
in terms of the error term £ijk as a scaled mixture of normal random variables 
to account for outlying observations. The error variance is defined in terms of 
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the measurement error variance erf., specific to cytotoxicity parameter j, and 
on ENM-specific variance inflation parameter n. If we define the joint distri- 
bution of £ ijk (d,t) and n as P(£ij k (d,t),Ti) = P(£ ijk (d,t) \ Ti,a Ej )P(ri \ v), 
choosing Eijk(d,t) \ Ti,a £j ~ N(0,a £ ./ri) and Tj | v ~ Gamma(^/2, z^/2), it 
can be shown that the marginal density of Eijk(d,t) \ a\. is distributed as 
a T((7j,,z/) [West (1984)]. Under this framework, we can borrow strength 
across all ENMs by assuming the error variance is the same, but retain ro- 
bustness in the model by allowing ENM-specific departures from normality. 
We allow the measurement error a £j to vary between cytotoxicity parameters 
due to heterogeneity in the cytotoxicity outcomes. 

Stage 2: Response model at the ENM by cytotoxicity parameter level. The 
dose-response surface rriij(d,t) spans two dimensions (dose and time), and is 
modeled in an additive fashion as described by Hastie and Tibshirani (1986). 
If we let (aij , /3'ij , 4>ij , 7y , iPij>d'ijiX'ij)' be a parameter vector indexing the 
dose-response surface mij(d,t), we can then define 

(3) mij(d,t) =a.ij + fijtd-^tjiPij) + gijit-^ij^ij) + hij(d,t;Xij,$ij)- 

Here fij(d;cj> i: j, (3^) is a function modeling the effect of dose d on response 
j for ENM i. Similarly, gij{t; 7/^,7^) is the function modeling the effect of 
time t and hij(d,t;Xij^ij) is the function modeling the interactive effect of 
dose and time. More specifically, we model the interaction of dose and time 
in a semi-parametric fashion as hij(dt; Xij^ij)- This parameterization allows 
us to retain direct interpretation of the model parameters, while avoiding 
over-fitting of sparse data. To ensure likelihood identifiability, we require, 
without loss of generality, that fij (d = 0; (j)^ , (3^ ) = 0, gij (t = 0; t/^.- , 7^ ) = 0, 
and hij(dt = 0; Xij,$ij) = 0- The parameters aij can therefore be interpreted 
as the background response level for each particle and outcome. 

We model dose-response curves fij (d; i • , ) , duration-response curves 
gij(t; tpij, 7^) and dose-time response curves hij{dt;Xij^ij) as linear com- 
binations of basis functions. Specifically, we use linear B-splines with two 
random interior knots as points where the slope changes in a piecewise lin- 
ear fashion. Let B(x, 77) denote a 4-dimensional B-spline basis with inte- 
rior knots 77 = (7/1,772)'. Also, let (3^ = (/%i,. . . ,/%4)', 7tj = (7ij'i> • • ■ Hiji)' 
and dij = (<%i, . . . , 5ij4)' be 4-dimensional vectors of spline coefficients. The 
functions f ij (d](f} ij ,f3 ij ), g^ (t; V^Tij) and hij{dt;Xij^ij) can then be rep- 
resented as follows: 

fij {d; 4>ij , /3 i:j ) = B(d, 4)^ YPij , 

(4) 9ij(t;^ij,lij) = S(t,Vii)%-, 
hij{dt; Xij^ij) = B{dt, XijY^ij- 

Identifiability restrictions fij(d = 0; 0^-,/3 ij -) = 0, gij{t = 0; 7/^,7^) = and 
hij(dt = 0;Xiji&ij) = are implemented by fixing = 0, 7^1 = and 
= 0, for all particles and outcomes (see Figure 2 for an illustration). 
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Fig. 2. Dose-response as a change-point model, (left) B-spline basis function of degree 1, 
corresponding to change points (interior knots) at log doses of 1.5 and 4-5. (Middle) Ex- 
ample dose-response curve. The basis function on the left corresponds to a spline function 
with 2 change points. Each random change point has a corresponding distribution, result- 
ing in a smooth dose-response curve. (Right) Example of a marginal prior distribution on 
the change points corresponding to the dose-response curve on the left. This formulation 
favors (a priori) the choice of conservative values for the location of the first change point 
(solid line), and a relatively diffuse prior for our second change point (dotted line). 

Modeling dose and duration-response curves as piecewise linear functions 
allows for considerable flexibility while maintaining direct interpretability of 
the model parameters. Recall that in our formulation the interior knots are 
estimated as random quantities. This allows, marginally, for a smooth dose- 
response trajectory that is automatically adjusted to fit the data. The main 
advantage of the proposed functional representation is that, in the absence 
of a dose-time interaction, one can interpret the first interior knot 4>iji as the 
dose at which ENM i becomes toxic in relation to cytotoxicity parameter 
j (Maximal Safe Dose — similar to the classical NOAEL concept). A similar 
interpretation can be given to ipiji, in relation to duration-response. Note 
that the foregoing interpretation is contingent on fixing = 0, 7^2 = and 
Xij2 = when assuming no effect before the first change-point, and /3y2 < 0, 
7ij2 < and Xij2 < when assuming a tonic effect before the first change 
point. In the presence of a dose-time interaction, interpretation changes 
slightly and we instead consider the idea of safe exposure regions, which 
represent doses and time exposure combinations that do not induce cyto- 
toxicity. Finally, in the absence of an interaction, the parameters ^>ij2 an d 
ipij2 are respectively interpreted as the dose and time at which the response 
stabilizes, or cells start a possible recovery process. 

We can expand the model further to allow for the exclusion of interaction 
functions where not needed. To do that, we include a latent indicator variable 
Pij, so that for each particle i and outcome j 

{aij + fij (d; 4>ij , f3 i:j )+gij (t; ^ , 7^ ) , if pij = 0, 
otij + fij(d; (hij-Pij) ■ H.jd-thj-lij) + hij(dt;Xij,8ij), 
if p i: j = 1, 
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where pij ~Bern(7r) and 7r~ [7(0,1). We require that if pij = 0, hij (dt; Xij , 
Sij) > 0, to ensure identifiability. The indicator variable p^ can then be used 
to test explicitly for the dose-time interactions. The exchangeable Bernoulli 
trials prior on p^j is designed to account for multiplicities [Scott and Berger 
(2006)]. This trans-dimensional parameterization is key to avoid overfitting, 
to facilitate parameter interpretation, and to allow for testing of specific 
scientific hypotheses related to the biological interference of nanomaterials. 

For each ENM i and response j, we define the following prior distributions 
for aij, (3^, j^, and 

aij~N(a 0i ,al.), 

(6) ^~Ar 4 (/3 o .,E^JJ{Aii = 0;/3y 2 <0;(^3,A^)>0}, 
-fij ~ JV 4 (7 .,£ 7t )I{7yi = 0;7i j2 < 0; (7^3, 7ij4) > 0}, 

Sij I pij = 1 ~ N^xaSi^vs^IiSiji = 0;<% 2 < 0; (<%3,<%4) > 0}. 

The truncated support for (3^, 7^ and Sij imposes functional constraints 
on /(■), g(-) and h(-), which are consistent with the expected behavior of 
canonical dose and duration kinetics. At the same time, however, it allows 
for the system to recover by permitting a decreasing slope after the second 
change point. The covariance matrix S^. has diagonal elements <rp , £ = 
1, . . . , 4, and off diagonal elements equal to 0; similarly for S~. . 

Prior distributions for (j)^ , and Xij are defined to satisfy the following 
constraints: (0 < fcji < 4>ij2 < D), (0 < ipiji < ^2 < T) and (0 < Xijl < 
Xij2 < DT). More precisely, we assume that the joint distribution of the 
interior dose and duration knots follows a generalized bivariate Beta density 
function, so that 

4>ij ~ B 2 (a^ 1 ,b ( j )1 ,a ( j )2 ,b ( j >2 ,D), 

(7) tpij ~ B 2 (a^ 1 ,b^ 1 ,a^, 2 ,b^ 2 ,T), 

Xij ~ B 2 (o xl , b xi , a X2 , b X2 , DT) . 

Here we assume that a random vector x = (sci, a; 2 )' is distributed according to 
a generalized bivariate Beta distribution function (x ~ B 2 {a\, b\, a 2 , b 2 , ni)), 
with support <S(x) = {(xi,x 2 ) : < x\ < x 2 < m} if and only if 

p(x I 01,61,02,62)"*) 

(8) =p(xi I ai,b 1 ,m)p(x 2 \ xi,a 2 ,b 2 ,m) 

1 xl 1 ~ 1 (m-xi) bl ~ 1 1 
~ B(ai,h) m ai+bi-i B(a 2 ,b 2 ) (m - x 1 ) a2+b2 ~ 1 

The foregoing formulation can be seen as a generalization of the Dirichlet 
distribution over a two-dimensional simplex. This general formulation can be 
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simplified further, in order to achieve a right-skewed marginal distribution 
for x\ and a uniform conditional distribution for x 2 given x\ . This is achieved 
by assuming b\ > a\ > 1 and a% = b 2 = 1. 

Making use of this construction, we simplify the prior distribution in (7) 
as follows: 

<kj ^-BaCljAfc^A^j,, 1, 1,D)I{ A^j, > > 1}, 

(9) ipij ~ B 2 (l,\f il ,\<p i2 ,l,l,T)I{\<i )i2 >)y Hl >1}, 

Xij ~ B 2 (l,l Xa ,l Xa ,lA,T)I{hca >l Xa > !}■ 
From a regulatory standpoint, this formulation favors (a priori) the choice of 
conservative values for the location of the first change point and a relatively 
diffuse prior distribution for our second change point (see Figure 2). 

Stage 3: Response model at the ENM level. For each ENM i, we exploit 
conditional conjugacy to define the following prior distributions for popula- 
tion level parameters: 

(10) a Qi ~N(m at ,v ai ), (3 . ~ N 4 (m ft , v ft ) , j . ~ N 4 (m^., v T .). 

In the absence of an interaction, the parameters /3 G . and *y . represent sum- 
maries of the dose and duration-response trajectories across all outcomes 
and the a Qi parameters represent a summary of the baseline response across 
all outcomes. In the presence of an interaction, we may construct these sum- 
maries conditionally on specific doses and durations of exposure. 

Finally, considering the distribution introduced in (9), we define a prior 
model for population level parameters A<£. = (A^ , A</, i2 ) and A^. = (A^ a , A^, i2 ) 
as follows: 

(11) A fe ~ Gamma(a A<M ,& v J, X^ u ~ Gamma(a^,6 A#( ), 

where £ = 1,2. The parameters A<^. and A^,. can be used to construct sum- 
maries of dose and duration-response change points across all outcomes. 
Shape hyperparameters , ) and (oa^A^J can be tuned to favor 
more or less conservative values for the change-point locations at the particle 
level. 

Stage 4: Hyperpriors. We complete the model by specifying prior distri- 
butions on our hyperparameters as follows: 

1/of ~ Gamma(o £j , b ej ) , l/o^. ~ Gamma(a a , , b a J, 

(12) 

1 / cr|. ~ Gamma(a / g i , b^) , 1/ <r^. ~ Gamma(a 7i , & 7i ). 

We model our precision parameters as gamma distributions, exploiting con- 
ditional conjugacy. Again, prior parameters can be tuned to define more 
or less informative distributions consistent with the scale of the outcomes 
[Gelman (2006)]. Note that in our formulation, x ~ Gamma(a,6) denotes a 
Gamma distributed random quantity with shape a and rate b, such that 
E(x) = a/b. 
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3. Estimation and inference. 



3.1. Posterior simulation via MCMC. Using the B-spline representation 
introduced in Section 2.1, we can write the expected jth response level 
associated with ENM i, at dose d and exposure time t as 



n 



+ B{d, ^YPij + B(t, %)'-Yij, if Pij = 0, 



mi(d, t; Oiij^ij,. ■■) = { OLij + B{d, fa^'fly + B(t, ip^'lij + B(dt, Xij)'8ij, 

if = 1. 

Let /3 = {/3^ : i = 1, . . . , /, j = 1, . . . , J} and define 7 and 5 in a similar 
fashion. These parameters denote the full set of spline coefficients. Further- 
more, consider knot parameters cp = {4>ij '■ i = 1> • • • > !■> j = 1) • • • 1 J}-> with tp 
and x similarly defined, and background response parameters o; = {ctij : i = 
1, . . . , /, j = 1, . . . , J}. Finally, let a 2 £ = {a 2 £l a^)' and r = (n, . . . , r/)'. If 
we denote with Y the complete set of response values for all particles and 
cytotoxicity outcomes, the likelihood function can be written as follows: 

L{fi,^,8,(j),ip,x,oc,crl,T,p\ Y) 

(13) 

'o- £ .\ I [ (y ijk (d,t) -m,ij(d,t;...)y 



« n 

i,j,k,d,t 



T i ) 6XP l 2 <r| /rj 



where the product is taken over all replicates k, particles i, outcomes j, 
doses d and times t. We are interested in the posterior distribution 

P(0,j,8,(p,ip,x,oc,(T%T,p\ Y) ocL(P,y,8,(p,ip,x,ac,(T 2 E ,T,p\ Y) 

(14) 

x P{/3,7,5,<j>,ip,x,ac,cr e ,T,p), 

where the prior model P(f3,~f,d,cj>,'ip,x, a i cr ' 2 i T i P) 1S IUU Y described in 
Section 2.1. This quantity is, however, unavailable in closed analytic form, 
therefore, we base our inference on Markov Chain Monte Carlo (MCMC) 
simulations. 

The proposed posterior simulation algorithm combines Gibbs steps within 
Metropolis-Hastings steps in a hybrid sampler, where we update parame- 
ters component- wise [Tierney (1994)]. We directly sample components when 
closed-form full conditional distributions are available using a Gibbs sam- 
pling algorithm [Geman and Geman (1984), Gelfand and Smith (1990)]; 
otherwise, we use the Metropolis-Hastings (MH) approach [Metropolis et 
al. (1953)]. Available full conditional distributions are given in the supple- 
mental article, Appendix A [Patel et al. (2012)]. As we are considering se- 
lection of interaction functions in a trans-dimensional setting, we implement 
a reversible jumps algorithm to move between models with and without the 
dose-time interaction function hij(dt;Xij,8ij) [Green (1995)]. The model in- 
dicator pij and corresponding model parameters Sij and Xij are updated 
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jointly using reversible jump MCMC steps. After the model structure has 
been specified, the model parameters are updated from their corresponding 
conditional posterior distributions. The proposed sampling scheme can be 
summarized as follows. 

1. Fixed dimensional updates. Given the current state of the latent in- 
teraction indicators py, response surfaces are uniquely defined as in (5). 
Posterior sampling is standard here and proceeds by updating spline co- 
efficients /3,7 and S from their conditional posterior via direct simulation 
[Patel et al. (2012)]. Knot parameters (f),if> and x are updated via a MH 
step. For example, when sampling the interior knot parameters (f) we use 
an appropriate proposal kernel q((/>%g, 4>\je) to efficiently construct Markov 
chains with the desired stationary distribution. While accounting for the 
fact that (fiiji < (f>ij2, we consider uniform proposal densities of the form 



where £ = 1,2. Here denotes the appropriate support and must satisfy the 
constraints < (j>ij\ < <pij2 < D. Proposed values of <f>iji are accepted with 
the following probabilities: 



To tune proposal kernels, each (ftiji was sampled using an initial value of w 
that was re-calibrated throughout the burn-in period to achieve an accep- 
tance rate between 30% and 70% [Roberts and Rosenthal (2001)]. Specifi- 
cally, the acceptance rate of (ftijg was monitored every 200 iterations through- 
out the burn-in period with adjusted appropriately if the acceptance 
rate did not fall within the desired range. A similar Metropolis-Hastings 
scheme was adapted for sampling the duration-response parameters dose- 
time interaction parameters Xij I Pij = 1, as well as for population level knot 
parameters. 

2. Trans- dimensional updates. We sample the model space by randomly 
proposing the birth or death of dose-time interaction functions hij(-). This 
is accomplished by selecting a particle i and outcome j at random and by 
jointly updating p^, Sij and Xij- I n detail, 

(1) For uniformly random i £ (1, . . . , /) and j S (1, . . . , J), propose a sys- 
tematic change p°- — > pjj = 1 — p?- . We assume for the moment that we pro- 
pose moving from p^ = to pjj = 1, implying the birth of a new interaction 
function /ijj(-). 

(2) Propose new knots and spline coefficients 8\j ~ q{8\j) and xh ~ 



(15) 



q(4>iji I 4>ijd = U (<p ije - w^ije, <f) ijt + w^jtjl^tf) 



(16) 




?(*«)■ 
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Table 1 

Risk assessment parameters. ENM level risk assessment parameters associated with the 
hierarchical model introduced in 2.1. For each parameter we summarize its function in 
the model and the related interpretation as a cytotoxicity risk factor 



Parameter 


Model function 


Parameter interpretation 




Dose-response slope from 4>tji to </>ij2 


Overall dose effect 




Duration-response slope from (j>iji to <\>iji 


Overall exposure time effect 


4>lij 


Dose-response change point 1 


Maximal safe dose 


Iplij 


Duration-response change point 1 


Maximal safe exposure time 




Evaluated numerically 


Maximal response 



(3) Accept the proposed move with probability r& = min(l,i?fc), where 

vivijk I p%,9\5 ihXij , Pij ) q&Mxjj) p(p°ij) ' 

where we use 6\ u to denote all model parameters, with the exception of oj. 

In the case where the proposed move would imply a death of an interaction 
function (p^ = 1 — > pjj = 0), the acceptance probability would simply be 

Td = l/n- 

While the proposal densities q(Sij) q(Xij) m (17) can in theory be defined 
almost arbitrarily, to guarantee efficient exploration of the model space, we 
consider truncated multivariate normal proposals for Sij and Xij centered 
around regions of high posterior probability. Efficient optimization within 
the MCMC iterations is achieved using standard profile likelihood ideas 
[Severini and Staniswalis (1994)]. 

3.2. Posterior inference. In this section we discuss inference on ENM- 
specific risk assessment parameters, based on draws from the posterior dis- 
tribution described in Section 3.1. Table 1 summarizes several quantities of 
interest including the maximal safe dose, maximal safe exposure time and 
the maximal response. This list is not exhaustive. However, other risk assess- 
ment parameters of interest, such as benchmark doses (BMD) or effective 
concentrations (ECa), are easily obtained from our model output in a nu- 
merical fashion. In the case of a dose-time interaction, these quantities are 
defined conditionally on specific doses and durations of exposure. 

T„ t An) ;,(n) (n) n {n) An) An) (n) , (n) _ -, AT A a 

Let (p^ , iptj , Xij , Pij , lij , , and p {j , n- l,...,iV, de- 
note N MCMC draws from the posterior distribution of (p^, tp^, Xij-, Piji 
ctij and p^. In the absence of an interaction term, posterior samples 

(n) (n) 

<Plj{ and tpljl directly provide us with an approximation of the posterior 
distribution for the maximal safe dose and maximal safe exposure time. 
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We can also obtain the posterior samples for the overall dose effect, = 

Piji/i^ijl ~ 4>ij\ )i which is the slope of the dose-response curve between 
and 4>ij2- Similarly, we can obtain the posterior distribution for the overall 

time effect using posterior samples 7^3 = 7*^3/(^2 ~' l Piji)- I* 1 the pres- 
ence of a dose-time interaction, we can define any of the summaries described 
above conditionally on a given dose and time. For example, the maximal safe 
dose conditional on exposure time can be defined as mm{(friji,Xiji/t}, and 

posterior samples can be obtained from min{0^,x^/t}. Given posterior 
draws, one can proceed with the straightforward construction of standard 
posterior summaries, such as means, maxima a posteriori, modes, quantiles 
and credible regions. We may also be interested in testing for a dose-time 
interaction. The expected inclusion probability of the dose-time interaction 
function can be estimated using posterior draws as pij = ^2 n p\^/N. 
Given the prior distribution described in (5), this posterior probability is 
known to adjust for multiplicities and can be used to test for a dose-time 
interaction. Scott and Berger (2006), for example, recommend selecting the 
median model, that is, including all interactions for which pij > 0.5. Also of 
interest is an estimate of the dose-response surface, rriij(d,t), for particle i 
and outcome j. This surface is, of course, defined in an infinite-dimensional 
space. However, given the basis-function representation introduced in Sec- 
tion 2.1, we only need finite draws from the parameter set of interest. More 
precisely, draws from the marginal posterior distribution of the dose-response 
surface for any dose d € [0, D] and time t € [0, T] are given by 



(18) m%\d,t) 



a£ n) + B(d, + B(t, » >, if p£> = 



\f +B( d ,^y^ +B(tM?)'j ( i : ) +B(dt,x% ) )'4?, 



■ r In) , 

ifp.y = 1. 



For each tjjjjy , ip^ , (3^ , <f>^ and offi , n = 1, . . . , N, we evaluate the dose- 
response function given in (18) over a grid of values D = (d\, . . . ,d n )' and 
T = (ti,...,t n )'. The posterior mean of the samples rn^\ n = 1,...,N, 

at each value of D and T can be used to summarize the fit of the dose- 
response surface, as shown in Figures 3 to 4. Other quantities of interest in- 
clude the posterior distribution of the dose-response function fij(d; 4> { j, (3^), 
duration-response function gij(t; ^p'yy) and dose-time interaction function 
hij (dt; Xij , Sij). Draws from the marginal posterior distribution of these func- 
tions for any dose d £ [0,D] and time t £ [0, T] are given by 

f^(d;<l> ij ,(3 ij )=B(d,^ ) )'^f, 

(19) 9^\t^l ij )=B{t^)'^\ 
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Fig. 3. Fitted response curves for the quantum dot ( QD) ENM. (left) Fitted response sur- 
faces (column 1), dose-response function, fij(d) (column 2), duration-response function, 
gij(t) (column 3), dose/duration interaction function, hij(dt) (column 4) and associated 
95% posterior intervals. In ( column 1 ), the color red represents response values corre- 
sponding to lower time points and the color black represents response values corresponding 
to higher time points. 



For each draw, we evaluate the dose-response functions over a grid of values 
d G D and the duration-response functions over a grid of values t £ T. As 
described before, standard pointwise posterior summaries can be obtained in 
a straightforward fashion. Simultaneous confidence bands for the functional 
effect of interest can be constructed following the Monte Carlo approxima- 
tion suggested by Baladandayuthapani, Mallick and Carroll (2005). 
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Fig. 4. Fitted response curves for the gold (Au) ENM. Fitted response surfaces (col- 
umn 1), dose-response function, fij(d) (column 2), duration-response function, gij(t) (col- 
umn 3), dose/duration interaction function, hij(dt) (column 4) o/nd associated 95% pos- 
terior intervals. In ( column 1 ), the color red represents response values corresponding to 
lower time points and the color black represents response values corresponding to higher 
time points. 



Additional summaries of interest can be obtained in a numerical fashion. 
For example, the posterior distribution for the maximal response value m*j = 

maxjmjj (d,t);d G [0, D],t E [0,T]}, may be obtained evaluating m\j\d,t) 
over a fine grid of doses D and times T. An approximate posterior draw from 
m*j can be defined as m*^ = max{roj" \d,t);d G D,t G T}. Given smooth- 
ness constraints on rriij(d,t), defined in Section 2.1, the foregoing procedure 
is likely to provide a good approximation to the posterior distribution of 
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the maximal response value, provided D and T define a sufficiently detailed 
evaluation grid. Similar procedures may be adopted to obtain inference on 
other risk assessment parameters like ECas or BMDs. 

4. Applications. 

4.1. Synthetic data. To assess estimation of the model presented in Sec- 
tion 2, we present a simulation study in the supplemental article, Appendix 
B [Patel et al. (2012)]. The dose and time kinetics were simulated in an 
additive fashion, from various parametric functions, including both canoni- 
cal and noncanonical profiles that are still reasonably interpretable under a 
toxicity framework. We also placed increasingly conservative priors on the 
population level parameters A^. and A<^. in order to assess the sensitivity 
of the model results to our choice of prior parameters. In the supplemental 
article, Appendix C, we provide an additional sensitivity analysis assessing 
model results to our choice of prior model for the change-point parameters 
[Patel et al. (2012)]. We compare our prior model results to both a truncated 
normal prior and a parameterization of the bivariate beta prior that results 
in a uniform prior on the simplex. 

Simulation results indicate that our model is robust to model misspec- 
ification and is not very sensitive to our choice of prior. We do, however, 
maintain that using the bivariate beta prior defined in (8) is likely to be 
more appropriate in data analytic frameworks, as the implied stochastic 
behavior of the response surface, a priori, reflects more closely the usual bi- 
ological mechanisms of toxicity. More specifically, it assigns zero probability 
of toxicity to zero dose and time, where toxicity is indeed not expected to 
occur. Furthermore, this prior accounts for issues such as dosimetry, in which 
the administered doses are confounded by different particle bioavailability. 
Therefore, in some particles toxicity is not expected to occur for doses and 
times greater than zero. 

4.2. Case study background. We illustrate the proposed methodology by 
analyzing data on macrophage cells (RAW cells) exposed to eight different 
metal and metal-oxide nanoparticles, monitored in relation to four cytotoxi- 
city parameters. All four outcomes are measured over a grid of ten doses and 
seven times (hours) of exposure (see Figures 3 to 4) . Cytotoxicity screening is 
based on the hierarchical oxidative stress model [George et al. (2009)]. More 
specifically, a multi-parametric assay that utilizes four compatible dye com- 
binations and the subsequent change in fluorescence read-outs was used to 
measure four responses relating to the highest tier of oxidative stress (toxic 
oxidative stress). The four measured responses include mitochondrial super- 
oxide formation (MSF), loss of mitochondrial membrane potential (MMP), 
elevated intracellular calcium (EIC) and cellular membrane damage (CMD). 
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Figure 1 provides fluorescence images of cells exposed to various nanoma- 
terials (50 fig/ ml and 3 hours), including quantum dot, platinum and a 
negative control consisting of no nanomaterials. Row 1 includes images of 
cells treated with a dye combination including MitoSox, which permeates 
the mitochondria and fluoresces red when oxidized by superoxide. Red flu- 
orescence measured in cells treated with MitoSox is therefore a measure of 
mitochondrial superoxide formation. Similarly, in Row 2 cells are treated 
with a dye combination including JC1, which stains the cytoplasm red in 
healthy cells, but forms a monomer in cells with decreased membrane po- 
tential and consequently stains the cytoplasm green. Finally, in Row 3 cells 
are strained with a dye combination including Fluo-4 and Propidium Iodide 
(PI). In cells with damaged membranes, PI is able to permeate the cell and 
bind to DNA where it causes the nucleus to emit a red florescence. Fluo-4 is 
a dye that emits a green fluorescence in the cytoplasm in cells with elevated 
intracellular calcium. Each sample was also stained with a Hoechst dye which 
causes all cell nuclei to emit a blue florescence, allowing for a count of the 
total number of cells. An analysis of the fluorescence readout, monitored at 
varying wavelengths, results in a measure of the percentage of cells positive 
for each response. Figure 1 also provides a heat map of the raw responses 
for each particle and outcome, where colder colors (blues and greens) in- 
dicate a smaller percentages of cells positive for the response and warmer 
colors (oranges and reds) indicate a higher percentage of cells positive for 
the response. The final data was normalized using a logit transformation to 
unconstrain the support so that it can take on values between — oo and oo. 
Our inferences are based on 20,000 MCMC samples from the posterior distri- 
bution in (14), after discarding a conservative 60,000 iterations for burn-in. 
MCMC sampling was performed in R version 2.10.0, and convergence diag- 
nostics were performed using the package CODA (Convergence Diagnostics 
and Output Analysis), [Plummerm et al. (2006)]. 

4.3. Case study analysis and results. We fit the model described in Sec- 
tion 2.1 to the metal-oxide data set described in the previous section. The 
prior on the interior knot parameters was modeled using the simplified den- 
sity described in (9). A set of relatively noninformative Gamma(2, 1) and 
Gamma(3, 1) priors were considered for the components of both A^. and 
A.0. , along with a vague i?2(2,3, 1, 1,DT) prior for our dose-time interaction 
change-point parameter Xij- We also fixed /3ij2 = and 7^2 = 0, assuming no 
effect before 0yi and ipiji, thereby allowing, in the absence of a dose-time in- 
teraction, the interpretation of <f>ij\ as the maximal safe dose and ipiji as the 
maximal safe exposure time. Similarly, when pij = 1, we fixed 6ij2 = 0. We 
placed Gamma(0.01, 0.01) priors on the l/c ej parameters, Gamma(l,0.1) 
priors on all remaining precision parameters and N(0, 10) priors on the a D i 
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Fig. 5. Graphical model diagnostics. (Bottom) Probability Integral Transform assessing 
empirical calibration of the posterior predictive distribution. (Top) Mean and 95% posterior- 
intervals of the posterior predictive mean response across all doses and times of exposure, 
for all outcomes and particles 1 through 8 (QD, ZnO, Fe^Oi, Pt, Ag, SiC>2, AI2O3, Au). 
Also included are the empirical mean responses across all doses and times of exposure 
(red). 



parameters. The parameters f3 oi and ~y oi are modeled as truncated multivari- 
ate normals with mean 1 and a covariance matrix with diagonal elements 
10 and off-diagonal elements 0. Finally, we placed a prior distribution on 
the degrees of freedom parameter v, for the T-distributed error described 
in Section 2.1. We specified the prior to be uniform on 1, 2, 4, 8, 16 and 
32 degrees of freedom [Besag and Higdon (1999)]. In concordance with our 
synthetic data experiments, a sensitivity analysis on the case study data set 
proved robust to reasonable variations in the prior specification. 

We provide graphical summaries of goodness of fit and posterior predic- 
tive performance in Figure 5. The top panel shows the mean and 95% poste- 
rior intervals of the posterior predictive mean response across all doses and 
times of exposure (black) , along with the empirical mean response (red) , for 
each particle and outcome. In all cases the empirical mean response is con- 
tained within the 95% posterior intervals of the posterior predictive mean 
distribution, indicating good average posterior coverage across doses and 
times of exposure. The bottom panel provides a plot of the probability in- 
tegral transform histogram for the entire model [Gneiting, Balabdaoui and 
Raftery (2007)]. Visual assessment of the plot indicates that it is close to 
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Fig. 6. Safe exposure regions for the quantum dot (QD) and platinum (Pt) nanomate- 
rials. For each particle and outcome we can define dose and time exposure regions which 
do not induce cytotoxicity. Red colored regions indicate greater cytotoxicity to the cells, 
whereas blue colored regions indicate reduced risk. Contour lines quantitate the median es- 
timated response, relative to the background, where zero response areas can be interpreted 
as safe exposure regions. 

uniformity, suggesting relatively good posterior predictive calibration. Ad- 
ditional summaries and diagnostic tools are detailed in the supplemental 
article, Appendix E [Patel et al. (2012)]. 

Figures 3 and 4 illustrate data and results associated with two of the 
particles examined in this HTS study. Particularly, we report inference for 
platinum and quantum dot nanomaterials for each of the 4 cytotoxicity 
outcomes. Inference for the remaining 6 particles is reported in the sup- 
plemental article, Appendix D [Patel et al. (2012)]. In these two figures, 
column 1 shows expected posterior dose-response surfaces across dose and 
time for all outcomes. As the posterior expectation marginalizes over the 
interior knots, smooth surfaces reflect the uncertainty about the location of 
these change points and provide an illustration of how the proposed tech- 
nique will adjust for smoothness in an unsupervised fashion. Also included 
are functional posterior expectations associated with dose-response curves 
fij(d) (column 2), which represent the effect due to dose, duration response 
curves gij(t) (column 3), which represent the effect due to exposure time, 
and the expected dose-time interaction function hij(t) (column 4). 

Figure 6 provides a plot of the estimated median response, relative to 
the background, for different doses and times of exposure. Blue colors in- 
dicate safety regions or areas of reduced risk to the cells, while red colored 
regions indicate increased risk of cytotoxicity. Finally, Figure 7 provides pos- 
terior summary estimates including mean and 95% posterior intervals for the 
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Fig. 7. Maximal Safe Dose for the quantum dot (QD) and platinum (Pt) nanomaterials. 
Posterior summary estimates of the maximal safe dose, conditional on exposure time, 
including the posterior mean and associated 95% posterior intervals. In the case of no 
interaction, the maximal safe dose is the same across all times. 



maximal safe dose, conditional on the duration of exposure. Note that in the 
absence of a dose-time interaction, the maximal safe dose is the same across 
all exposure times. 

Quantum dot (QD) shows a relatively high toxic response for plasma 
membrane damage and mitochondrial superoxide formation. In particular, 
we see a more pronounced dose effect for membrane damage and both a 
time, dose and significant dose-time interaction (p = 0.99) effect for mito- 
chondrial superoxide formation. This supports what has previously been 
demonstrated in conventional assays that QD nanoparticles stabilized by 
toluene are capable of inducing tiers 2 and 3 oxidative stress responses in- 
duced by the toluene [George et al. (2011)]. Platinum (Pt) shows a high 
dose and time response for mitochondrial superoxide formation, including 
a significant dose-time interaction effect (p = 0.99), and a pronounced time 
effect for elevated calcium but not for mitochondrial depolarization or mem- 
brane damage, indicating that the particle induced sublethal effects to the 
cell without cytotoxicity. The Zinc oxide nanoparticle (ZnO), reported in the 
supplemental article, Appendix D, shows a relatively high toxic response for 
plasma membrane damage, elevated calcium and mitochondrial depolariza- 
tion [Patel et al. (2012)]. In particular, we see a more pronounced time effect 
for the elevated calcium and both a time and dose-response for membrane 
damage and mitochondrial depolarization. This again verifies what has pre- 
viously been demonstrated in conventional assays, since ZnO nanoparticles 
are capable of inducing tiers 2 and 3 oxidative stress responses through Zn^~ 
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release [George et al. (2009)]. In contrast, the gold nanoparticle (Al), also 
reported in the supplemental article, Appendix D, shows very little response 
for all outcomes, indicating that, compared to the other particles, it has small 
risk of inducing a sublethal or lethal cytotoxic response [Patel et al. (2012)]. 

5. Discussion. In this article we propose a statistical framework for mod- 
eling dependent dose-response surfaces over multivariate outcomes. The pro- 
posed methodology accounts for dose and duration kinetics jointly using a 
flexible model which does not compromise interpretability. We account for 
the multivariate nature of the data using the hierarchical framework and 
thereby efficiently combine information and borrow strength across cellular 
injury patterns. We account for the nonrobust nature of the data by allowing 
for particle specific variance inflation, resulting in a T-distributed model for 
the error structure. 

The main challenge associated with the class of models proposed in this 
manuscript is finding the right balance between model complexity and model 
interpretability. An alternative formulation of the dose-response surface would 
seek inference for a general smooth surface rriij(d,t). However, our simplified 
approach, based on the assumptions of additivity and linearity, maintains a 
very appealing level of interpretability, allowing for the definition of specific 
risk assessment parameters while maintaining an adequate level of flexibility. 
A related generalization of the proposed additive framework would include 
a more general class of functional interactions to account for a possible syn- 
ergistic effect between dose and duration of exposure. This would come at 
the cost of reduced interpretability, but, at the same time, could be of clear 
scientific interest in some contexts. In this initial modeling effort, we choose 
to work with a T-distributed error structure and therefore normalize our re- 
sponse to unconstrain the support so that it can take on values between -co 
and oo. An alternative formulation could retain the original scale of the data, 
but rather define a generalized multivariate model such that the outcome 
distribution can be described using binomial or beta random quantities. This 
extension would require a substantial increase in computational complexity, 
with the possible need to consider numerical or analytical approximations, 
but it is clearly worthy of further methodological exploration. 

The hierarchical formulation introduced in this article is easily adapted to 
the case where multiple cell lines are used to test for cytotoxicity. A natural 
integration strategy would perhaps find motivation in the meta analytic 
framework, with information shared between experiments via the structuring 
of one extra level in the hierarchy. 

Finally, the proposed model can also be expanded by the inclusion of 
covariates. This is naturally defined as an extension to stage 3 of the model 
introduced in Section 2. The addition of covariates is especially important for 
relating specific ENM properties to toxicity, and is therefore an important 
area for future work. 
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SUPPLEMENTARY MATERIAL 

Supplementary Appendices (DOI: 10.1214/12-AOAS563SUPP; .pdf). Full 
conditional distributions for the model described in Section 2 are provided 
in the supplemental article, Appendix A. Spline coefficients /3,7 and S are 
directly sampled from their conditional posterior distributions via direct 
simulation (Gibbs step). To assess estimation of the model presented in 
Section 2, we present a simulation study in the supplemental article, Ap- 
pendix B. The dose and time kinetics were simulated from various paramet- 
ric functions. Both canonical and noncanonical profiles that are reasonably 
interpretable under a toxicity framework were generated. In addition, we 
assess sensitivity of the model results to our choice of prior parameters for 
population level interior knot parameters A<£. and A^ . In the supplemental 
article, Appendix C, we provide an additional sensitivity analysis assessing 
model results to our choice of prior model for the change-point parameters. 
Alternative prior models assessed include a truncated normal prior and a 
parameterization of the bivariate beta prior that results in a uniform prior 
on the simplex. The supplemental article, Appendix D, presents results as- 
sociated with inference on the 6 remaining particles not presented in Sec- 
tion 4.3. Finally, Appendix E discusses model assessment and goodness-of-fit 
diagnostics associated with the model described in Section 2. 
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